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Many soft-matter and biophysical systems are composed of monomers which reversibly assemble 
into rod-like aggregates. The aggregates can then order into liquid-crystal phases if the density is 
high enough, and liquid-crystal ordering promotes increased growth of aggregates. Systems that 
display coupled aggregation and liquid-crystal ordering include wormlike micelles, chromonic liquid 
crystals, DNA and RNA, and protein polymers and fibrils. Coarse-grained molecular models that 
capture key features of coupled aggregation and liquid-crystal ordering common to many different 
systems are lacking; in particular, the role of monomer aspect ratio and aggregate flexibility in 
controlling the phase behavior are not well understood. Here we study a minimal system of sticky 
cylinders using Monte Carlo simulations and analytic theory. Cylindrical monomers interact pri- 
marily by hard-core interactions but can stack and bind end to end. We present results for several 
different cylinder aspect ratios and a range of end-to-end binding energies. The phase diagrams are 
qualitatively similar to those of chromonic liquid crystals, with an isotropic-nematic-columnar triple 
point. The location of the triple point is sensitive to the monomer aspect ratio. We find that the 
aggregate persistence length varies with temperature in a way that is controlled by the interaction 
potential; this suggests that the form of the interaction potential affects the phase behavior of the 
system. Our analytic theory shows improvement compared to previous theory in quantitatively pre- 
dicting the I-N transition for relatively stiff aggregates, but requires a better treatment of aggregate 
flexibility. 
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I. INTRODUCTION 

Self-assembly of aggregates is ubiquitous in soft matter 
and biophysics [TJ. Aggregation requires only a pool of 
monomers with some type of attractive interaction. Re- 
versible aggregation occurs when monomers interact via 
relatively weak, non-covalent attractions that are typical 
of soft matter; in this case the aggregates are in equilib- 
rium with the pool of monomers. Reversible aggregation 
is also called equilibrium polymerization or supramolcc- 
ular polymerization [2; . While different geometries of ag- 
gregates are possible (depending on the nature of the at- 
traction) , in many cases anisotropic interactions favor lin- 
ear or filamentous aggregates. These rod-like aggregates 
can form liquid-crystal (LC) phases; the liquid- crystal 
order then couples to the aggregation, often promoting 
the formation of longer aggregates in the LC phase (H 0] . 
The liquid-crystal ordering of the aggregates can occur 
even when the monomers alone do not form liquid-crystal 
phases. 

Since the ingredients required for coupled aggregation 
and LC ordering are so simple, this basic physics oc- 
curs for a variety of systems, including worm- like mi- 
celles and microemulsions, chromonic liquid crystals, nu- 
cleic acids, proteins, and protein assemblies. Some of 
the first work on coupled aggregation and LC order was 
inspired by experiments on sickle-cell hemoglobin pro- 
tein, a fiber-forming protein important in sickle-cell ane- 
mia [5J [B] ; later work has considered other peptides and 



proteins that form fibers or fibrils [71H^]. Other early 
work considered worm-like micelles, formed either of am- 
phiphilic molecules in water or microemulsions of water 
and oil which are stabilized by amphiphilic molecules. If 
the micelles are relatively stiff they can display liquid- 
crystal phases [13]; these systems have been the subject 
of extensive experimental study [LlrES] . Liquid-crystal 
phases have also been observed in folic acid salts [29] . 
guanosine derivatives E01 EI] ; and nucleosome core par- 
ticles [32] . Recent work on chromonic liquid crystals and 
nucleic acids has renewed interest in coupled aggregation 
and LC order. Chromonic LCs are formed from relatively 
flat, disk-like dye molecules. When the dye molecules 
are suspended in water, hydrophobic interactions drive 
the stacking and formation of rod-like aggregates that 
can form LC phases [53TH5] . Hydrophobically-driven end 
stacking also causes short pieces of DNA and RNA to as- 
semble into aggregates and form LC phases [46H48] . 



Extensive work on the analytic theory of aggregation 
and LC ordering has been done since the early work of 
Herzfeld and Briehl [5 . Much of the initial work focused 
on perfectly rigid aggregates [SIlrJI HMMj . However, when 
aggregates are assumed perfectly rigid and excluded vol- 
ume interactions are treated in the second-virial approxi- 
mation, a "nematic catastrophe" occurs in which the ag- 
gregates in the nematic phase grow infinitely long [551I56J . 
Later work therefore emphasized adding aggregate semi- 
flexibility to the models [57H66] . Analytic theory qual- 
itatively reproduces results of experiments and Simula- 



tions, including a first-order I-N transition and longer 
aggregates in the nematic phase. Some work has also 
predicted an isotropic-nematic-columnar triple point |54j . 
However, quantitative agreement between simulation and 
theory remains lacking, particularly for aggregates with 
higher flexibility [65j [66] • 

Surprisingly little simulation work has addressed cou- 
pled aggregation and LC order. Some papers have 
focused on the molecular details that promote aggre- 
gate formation in specific systems and considered single 
aggregates [67- 72 . Simulating larger systems of interact- 
ing aggregates is too computationally costly for atom- 
istic simulation and coarse-grained models are required. 
Edwards, Henderson, and Pinning did pioneering Monte 
Carlo simulations of disks formed from hard spheres, with 
interactions that promoted disk stacking. They observed 
isotropic, nematic, columnar, and crystalline phases [73] . 
A similar approach was used by Maiti et al. to simu- 
late formation of columnar aggregates in chromonic liq- 
uid crystals. This work focused on conditions required 
to find columnar aggregates and didn't study the phase 
behavior [71]. Rouault used a coarse-grained 2D lattice 
model of wormlike micelles to study how varying the stiff- 
ness affects the ordering transition; stiffer molecules dis- 
play higher orientational ordering |75| . 

More recently, a number of simulation papers have 
considered hard spheres with an added anisotropic in- 
teraction that promotes linear aggregation. Hentschke 
and coworkers performed molecular dynamics and Monte 
Carlo simulations of spheres with anisotropic Lennard- 
Jones potential (76] [77]. They observed isotropic, ne- 
matic, and columnar phases, and observed the narrow- 
ing and eventual disappearance of the nematic phase in a 
I-N-C triple point as the temperature is increased. Simi- 
larly, Chatter ji and Pandit added an anisotropic interac- 
tion to hard spheres via a 3-particle potential that favors 
linear aggregation 78 . They observe a first-order I-N 
transition in Monte Carlo simulations. Some of the most 
detailed simulation work of a simplified model was done 
by Lii and Kindt, who performed Monte Carlo simulation 
of spheres with sticky "patches" that allow assembly into 
linear aggregates 65, 66. Their simulations did not show 
the columnar LC phase but did see the I-N transition. 

The key physical ingredients required for aggrega- 
tion coupled to LC order are (1) monomers with an 
anisotropic attractive interaction that promotes the for- 
mation of long, thin filaments and (2) excluded volume 
interactions between aggregates that promote liquid- 
crystalline order at high aggregate density. While other 
physical effects such as electrostatic interactions are 
clearly important in some systems [7J [T2] [251 E21 175] - 
a minimal physical picture which incorporates filament 
formation and excluded- volume interactions is a valuable 
starting point to understand common features of the di- 
verse experimental systems. 

In this work we develop a coarse-grained model of 
aggregation-induced LC order, and study the model us- 
ing Monte Carlo simulations and analytic theory. Our 



goal is to ultimately compare the simulation and the- 
oretical results to experimental data on chromonic and 
DNA liquid crystals, so we have chosen the form of the 
model and the approximations used to facilitate this com- 
parison. Motivated by recent work on nucleic acids and 
chromonic liquid crystals, we consider molecules with 
anisotropy both in shape and interactions. Therefore, 
our monomer is a cylinder of length L and diameter D, 
so the aspect ratio r = L/D. The cylinders are sticky: if 
two cylinders stack so that their circular ends are suffi- 
ciently close together, they experience an attractive inter- 
action. Other than this binding energy, the cylinders ex- 
perience only a hard-core repulsion. This coarse-grained 
model captures the main physical effects governing linear 
aggregation and mesophase formation in these systems, 
namely excluded volume interactions and end-to-end as- 
sociation. For the appropriate range of binding energy, 
density of cylindrical monomers, and temperature, the 
Monte Carlo simulations of sticky cylinders find linear 
aggregates and a range of phases, including isotropic, 
nematic, columnar liquid crystal, and columnar crystal 
phases. 

The hard cylinder system (corresponding to zero bind- 
ing energy between cylinders) has been studied previ- 
ously in simulations by Blaak et al. [50] and by Ibarra- 
Avalos et al. [81], but to our knowledge this is the first 
simulation study of hard cylinders with attractive inter- 
actions. 

In simple (non-aggregating) LC systems, the rod as- 
pect ratio is a primary determinant of possible phases 
and the location of the phase transitions. The role of 
monomer aspect ratio in controlling the phase behav- 
ior of aggregating LC systems is not well understood. 
Our sticky-cylinder model enables systematic investiga- 
tion of the dependence of aggregation and phase behavior 
on monomer aspect ratio L/D for systems ranging from 
disk-like chromonic liquid crystals (L/D -C 1) to duplex 
DNA oligomers (L/D > 1). 

In section [H] we discuss the simulation model and the 
Monte Carlo simulation methods, including cluster moves 
we used to improve configurational sampling. Section [III] 
outlines the analytic theory, the free energy minimiza- 
tion, and the I-N phase coexistence equations. In section 
|IV| we present the simulation results, including the phase 
behavior, order parameter and pair correlation functions, 
aggregate length distributions, and aggregate persistence 
length. The comparison of analytic theory and simula- 
tion for the I-N transition is presented in section [V| with 
a comparison both to our simulations and to the work 
of Lii and Kindt [S5]. Section VI is the conclusion. Ap- 



pendix [A] describes some calculation details of the ana- 
lytic theory. 



II. SIMULATION METHODS 

The monomers are hard cylinders of length L and di- 
ameter D with short-range attractive interactions be- 



tween cylinder ends. The short-range attraction between 
cylinder ends is a generalized ramp potential, 



U(r) 



0, 



-E hond [1 - {r/r c )~ 



r <r c 
r > r r . 



(1) 



where r is the distance between the centers of the circu- 
lar end faces of two neighboring cylinders. An exponent 
of 7 = 1 corresponds to a linear ramp potential, while 
the limit 7 — ¥ 00 corresponds to a square well potential. 
We chose 7 = 2 and r c — D/2, parameter values that 
we empirically found promote the formation of linear (as 
opposed to branched) aggregates; larger values of r c and 
7 give rise to branched networks. 

This model depends on three dimensionless parame- 
ters: cylinder aspect ratio L/D, packing fraction (f> = 
vq/v (or, alternatively, dimensionless pressure (3PD 3 ), 
and dimensionless binding energy /3£"bond- Here P is 
the pressure, /3 = l/(fcgT) is the inverse temperature, 
t>o = nD 2 L/A is the cylinder volume, and v = V/N is 
the volume per particle (inverse number density). Unless 
otherwise noted, the cylinder diameter D is the unit of 
length. 

We investigated systems of sticky cylinders with three 
aspect ratios, L/D = 0.5, 1, and 2. We carried out NPT 
Monte Carlo simulations over a range of pressures for 
binding energies ranging from /3-Ebond = to /?-Ebond = 
12. 

The persistence length l p of sticky cylinder aggregates 
is implicitly determined by the cylinder-cylinder pair in- 
teraction potential; no explicit bond angle bending po- 
tential is included. An implicit dependence of l p on the 
nature and strength of effective intermolecular interac- 
tions is a key feature of nucleic acid and chromonic ag- 
gregates, which may influence the temperature depen- 
dence of the I-N coexistence curve. In this respect our 
model is more specific than that of Lii and Kindt |55] . 
which included a bond angle bending potential, enabling 
the binding energy and persistence length to be varied 
independently. 



A. Monte Carlo Simulations 

We carried out Monte Carlo (MC) simulations of sticky 
cylinder systems with periodic boundary conditions, with 
fixed monomer number N, pressure P, and temperature 
T JS2J. Our simulations employed several types of MC 
move, including small displacements and rotations of sin- 
gle particles and changes in simulation cell volume and 
shape under an applied pressure. To improve configura- 
tional sampling, we also used cluster moves (see below), 
and flip moves (7r/2 reorientation of individual cylinders). 
An MC cycle includes N trial single-particle displace- 
ment/rotation moves, one trial volume-changing move, 
and (in cases where these are used) N/10 trial cluster 
moves and/or AT/10 trial flip moves. 

We simulated systems with between N — 720 and 
N = 1920 cylinders. For each aspect ratio and binding 



energy we performed both an expansion run (a series of 
simulations for decreasing pressure starting from a high- 
density columnar crystal) and a compression run (a series 
of simulations for increasing pressure starting from a low- 
density isotropic fluid). For each value of the pressure, 
we simulated 10 6 MC cycles, and measured thermody- 
namic and structural properties during the final 5 x 10 5 
MC cycles. To accommodate crystalline and LC phases 
of arbitrary symmetry without defects or strain, we uti- 
lized a fully flexible simulation cell, except in cases where 
this led to highly elongated or deformed simulation cells 
(e.g., at low densities). 



1. Cluster moves. 

Long aggregates undergo slow effective translational 
and rotational diffusion in simulations that utilize only 
single-particle MC moves. To improve configurational 
sampling of aggregates, we used cluster moves: simulta- 
neous displacements and rotations of groups of neighbor- 
ing cylinders. Similar cluster moves have been used in 
simulations of ionic fluids [55]. Clusters are defined by 
proximity: any two cylinders with r < r c i us t belong to 
the same cluster. The cluster cutoff r c i us t can be var- 
ied to control the average cluster size (n c iust) ■ For the 
simulations described here, we adjusted r c i ust to main- 
tain (n c i us t) ~ 3 for all pressures, binding energies, and 
aspect ratios. 

A trial cluster MC move consists of the following five 
steps. 1. Select a root cylinder at random, and find all 
cylinders in the cluster containing the root cylinder. 2. 
Generate a trial displacement and rotation of the cylin- 
ders in the cluster. The cluster is rotated as a rigid body 
about an axis passing through the center of the root cylin- 
der. We performed two types of rotation move with equal 
probability: spin moves, for which the axis of rotation is 
the axis of the root cylinder, and tilt moves, for which 
the axis of rotation is a randomly selected axis perpen- 
dicular to the root cylinder axis. 3. Reject the trial 
move if cylinder-cylinder overlaps are found in the new 
configuration. 4. Reject the trial move if the number 
of cylinders in the cluster has changed (i.e., if cylinders 
not in the original cluster are members of the cluster in 
the new configuration), as such moves violate detailed 
balance. 5. If no overlaps are found, and if the clus- 
ter size is unchanged, accept the move with probability 
V = min[l,exp(-/3Al/)], where AC/ = f7 fina i - initial is 
the change in potential energy. 

Cluster moves were used in all simulations of cylin- 
ders having a nonzero binding energy (/3-Ebond > 2). 
Flip moves were previously introduced by Blaak et al. 
[80] to improve configurational sampling for hard cylin- 
ders at high density, in the vicinity of a partially ordered 
(cubatic-like) intermediate phase. We used flip moves for 
small binding energy (-Ebond < 2) for the same reason. 

Cluster moves moderately improve sampling efficiency, 
as measured by the cylinder orientational decorrelation 
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FIG. 1. Orientational autocorrelation function C(t) for an 
isotropic fluid of sticky cylinders with aspect ratio L/D — 0.5 
and binding energy /3-Ebond = 12 at a pressure of f3PD 3 — 0.8. 
Results are shown for simulations with (solid curve) and with- 
out (dashed curve) cluster MC moves. Cluster moves decrease 
the orientational decorrelation time by approximately a factor 
of two. 



rate. The orientational autocorrelation function is 



C(r) 



1 N 
t=i 



[(iii(t) ■ Ui{t + t)) - (Ui(t)) • (u<(t + t))] , 

(2) 

where u^ is a unit vector along the axis of cylinder i, and 
the angle brackets denote an average over all time origins 
t (with "time" measured in units of MC cycles). In Fig- 
ure [I] we show the orientational autocorrelation function 
C(t) for an example simulation in the isotropic phase. In 
this simulation, cylinders assemble into aggregates con- 
taining about 10 cylinders on average, which leads to 
slow effective rotational diffusion. Note that C(t) de- 
cays exponentially for large r. Cluster moves reduce the 
orientational correlation time by approximately a factor 
of two, from 7.6 x 10 4 MC cycles in a simulation with only 
single-particle moves (dashed curve) to 3.8 x 10 4 MC cy- 
cles in a simulation incorporating cluster moves (solid 
curve). This decrease occurs despite the fact that the 
acceptance rate for cluster moves is rather low (around 
3.5% in this case). Performing TV/10 cluster moves per 
MC cycle adds approximately 50% to the computational 
cost of the simulation, so the gain in efficiency is around 
35%. Similar performance was found for other parame- 
ter values. More substantial improvements in sampling 
efficiency will likely require more sophisticated schemes 
such as configurational bias Monte Carlo [82] or the clus- 
ter cleaving method of Whitelam and Geissler [8H US] ■ 



2. Overlap test. 

The test for cylinder-cylinder overlaps is more time- 
consuming than the corresponding test for hard spheres, 
spherocylinders or ellipsoids. We use a procedure similar 



to that used in previous simulation studies of hard cylin- 
ders [80j El], and that consists of three tests: 1. Sphe- 
rocylinder overlap: test for overlaps between the enclos- 
ing spherocylinders (obtained by adding hemispherical 
endcaps to hard cylinders). If the spherocylinders don't 
overlap, then neither do the hard cylinders. 2. Disk-disk 
overlap: if an overlap is found in the first step, test for 
overlaps between the flat end faces of the cylinders. 3. 
Disk-cylinder overlap: if no overlaps are found in the sec- 
ond step, test for overlaps between the flat end faces of 
one cylinder and the cylindrical surface of the other. 

If the enclosing spherocylinders overlap, but no 
cylinder-cylinder overlaps are found, then the two cylin- 
ders may interact via the short-range attractive potential 
(Eq. (U])). Computing the attractive interaction between 
sticky cylinders requires minimal additional effort, be- 
cause the distance between the centers of the circular end 
faces of the cylinders is available from the spherocylindcr 
overlap test. 



III. THEORY 

To compare the I-N transition in simulation and the- 
ory, we extended the work of van der Schoot and Cates 
[56) . This analytic theory considers cylindrical monomers 
(of length L and diameter D) that reversibly assemble 
into linear aggregates. The form of the free energy is 
determined from three key assumptions. First, the ag- 
gregates are treated as nearly rigid cylinders. The ratio 
of aggregate length £ to the persistence length l p is as- 
sumed small, £/l p <C 1. Second, the aggregates interact 
via steric repulsion. Only pairwise interactions are taken 
into account within the second virial approximation, but 
the Parsons-Lee approximation is used to extend the the- 
ory to higher density. Third, the assembly of monomers 
into aggregates is driven by an energetic penalty -Ehond 
for monomers to have free ends. We assume that -Ebond 
is independent of the aggregate length £. 

The assumption of nearly rigid aggregates does not 
hold in our simulations: the results shown in this paper 
have a ratio of aggregate lengths to persistence length 
£/l p ~ 1. We discuss below the limitations of the theory 
in this flexibility regime. 

The work of van der Schoot and Cates considered two- 
particle interactions in the second virial approximation. 
Since this approximation is accurate only in the dilute 
limit, we used the Parsons-Lee approximation to extend 
the theory to the regime of higher packing fractions [5j5] . 
The Parsons-Lee method effectively inlcudes contribu- 
tions from higher virial coefficients in the interaction free 
energy. 

The total free energy of the system depends on the 
density of monomers pe(u) that belong to a subset of ag- 
gregates of contour lengths between £ and £+d£ and have 
orientations within the solid angle f2 u and f2 u + c?J7 u with 
respect to some specified direction n. The monomers 
have length L, diameter D, aspect ratio r = L/D, and 



volume vq = ttD 2 L/4:. The free energy per unit volume 
(in units of fcgT) is 
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Jdl— lPe (u) 
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4tt 



W(u). (3) 



All length integrals (over £) are from zero to infinity and 
all angular integrals (over f2) are over Air of solid angle. 
The first term is the free energy of a mixture of ideal 
gases of aggregates. Aggregates of different lengths and 
orientations are treated as different chemical species. The 
second term describes semiflexibility of aggregates; it is 
a perturbative correction to the orientational degrees of 
freedom (implicitly included in the first term) due to the 
flexibility of aggregates [57] . The third term describes the 
excluded- volume interaction between pairs of monomers 
with relative orientations described by cos 7 = ii! • u 2 . 
The prefactor i](<fi) reflects the Parsons-Lee correction 
for higher virial coefficients [861 US US] based on the 
Carnahan-Starling expression for the correlation function 
of hard spheres [90] 1 



T)(4>) 



1 / 4-34) 

4 {a^w 



(4) 



where <f> = pvQ is the packing fraction of monomers. In 
the usual second- virial approximation r](4>) = 1. The 
last term in Eq. (pi) is the free energy of polymerization, 
which is proportional to the total number of aggregates 
in the system. In writing Eq. ([3]) we assumed that in 
equilibrium the distribution of monomers is spatially ho- 
mogenous and that the conformational degrees of free- 
dom of aggregates decouple from the interaction degrees 
of freedom. 

For a closed system the total number of monomers N 
is fixed, leading to the normalization condition 



M ^P^) 



N 

V 



(5) 



To determine the I-N phase diagram we take the stan- 
dard approach of requiring mechanical and diffusive equi- 
librium of the two phases at coexistence. We calculate 
the osmotic pressure •p^ l ' n > and chemical potential per 
monomer p( l ' n > in each phase. Coexistence requires 



pd) = p (v) t 



n® = n<- n \ 



(6) 



The solution of Eqs. ffih provides information about the 
packing fractions <fi and <jy- n ' of monomers at the co- 
existence boundaries, the mean aggregation number (n) 
(and distribution of aggregation numbers), as well as the 



order parameter S = |((3cos 2 0) — 1), where is the an- 
gle between the monomer's axis and the director fi, and 
averaging is performed over all monomers in the system. 
In the sections below we outline the calculations and 
resulting equations for the isotropic and nematic phases 
as well as the numerical methods. 



A. Isotropic phase 

In the isotropic phase the monomers are distributed 
uniformly both in space and in orientation, so pf (u) = 

(i) 

p\ independent of angle. The free energy functional 
Eq. (pf for the total free energy of the system in the 
isotropic phase reduces to 
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where the monomer aspect ratio r = L/D and the 
monomer packing fraction = pvo. 

Functional minimization of the free energy is subject 
to the normalization condition Eq. (|5| , which can be han- 
dled with a Lagrange multiplier A: 
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This results in a length distribution of aggregates that is, 
as assumed, independent of orientation 
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(10) 



This solution corresponds to an exponential distribution 
for the density of aggregates Na gg (£)/V as a function of 
aggregate length t. 
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The chemical potential per monomer is given by 
» {l) = ^rr=rct>{2r,{<l))+H{<t>)) 



(11) 
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where r]'((f>) is the derivative of rj((j)) with respect to 
The osmotic pressure is 
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FIG. 2. Instantaneous configurations from NPT MC simulations of TV = 1440 L/D = 2 sticky cylinders with binding energy 
jfl-Ebond = 12 in various phases: (a) columnar LC phase, f3PD :i = 0.7; (b) nematic phase, f3PD z — 0.28; (c) isotropic fluid 
phase, PPD 3 = 0.18. Cylinder color encodes orientation: the RGB color index of cylinder i is (|«ia:|, |tliy|, |t(fa|), 
unit vector along the cylinder axis. 



where u; is a 
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FIG. 3. Equation of state for L/D — 1 measured on expansion from a high-density columnar crystal state (symbols). The 
horizontal lines indicate the approximate location of first-order phase transitions. 



B. Nematic phase 

In the nematic phase the system is spatially homoge- 
nous but orientationally ordered due to alignment of ag- 
gregates along a spontaneously chosen direction n. The 
system preserves azimuthal symmetry with respect to ro- 
tations about n. Therefore, in spherical coordinates the 
density pe{u) is independent of the polar angle <p and 



dO u , 1 

Pf[U) = - 

4tt m ' 2 



dxpt{x), 



(14) 



with x — cos 9, where 9 is the angle to the director n. 
As in the previous section, we minimize the system's free 
energy with respect to the orientation-dependent density 



function pi(x), 

SiFM+Xfdif^dxpiix)) 
5pt{x) 

This results in the integro-differential equation 
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The interaction kernel is 
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Solving Eq.(16) for the function pe(x) that minimizes 



the free energy is a difficult task that requires either nu- 
merical solution or introduction of a trial function. We 
chose a mixed approach that starts with a trial func- 
tion and uses numerical solution of the resulting alge- 
braic equations. We adopted the trial function proposed 
in reference 1561: 
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where the parameter Mq is given by Eq. ( 10 1 and M(x) 
is an orientation-dependent aggregation number. 

The free energy of the nematic phase as a functional 
of M(x) becomes 
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The bars denote dimcnsionlcss variables: l p — l p /L and 
M(x) = M(x)/Mq. The next step is functional min- 
imization of the free energy to determine an integro- 
differential equation for M(x). We find a series solution 
for M(x) by expanding in Legendre polynomials, which 
turns the equation for M(x) into a system of algebraic 
equations that we solve numerically. For details of the 
calculations, see Appendix A. Here we summarize our 
results for the chemical potential per monomer fjM 1 ' and 
osmotic pressure p^ in the nematic phase: 
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where I1—I3 are integrals involving M{x) that are defined 
inEqs. dA5|) - dA7 ) . 



The mean aggregation number in the nematic phase is 



(n) 



fdeN™w 



where 



Ng>W = (L/e) / pf^)^ 



(22) 



(23) 



is the distribution function of aggregates of length £ in 
the nematic phase. The order parameter of aggregates of 
length £ is 



St = 



§_ 1 dxP 2 {x)p l {x) 
f_j dxpt(x) 



(24) 



where ^2(2;) is the second-order Legendre polynomial. 
The order parameter (averaged over all monomers in the 
system) is 



J d£f_ 1 dxP 2 (x)pe(x) 
J d£ / 1 dxpi(x) 



(25) 



IV. SIMULATION RESULTS 

In Figure [2] we illustrate the phase behavior with sam- 
ple snapshots of the simulations. Instantaneous configu- 
rations from simulations of the L/ D = 2 sticky-cylinder 
system with binding energy /3-Ebond — 12 are shown. This 
figure displays the typical phase sequence observed upon 
expansion from a high-density columnar crystal for large 
binding energy, which includes the columnar LC phase, 
the nematic phase, and the isotropic fluid phase. 

To determine phase boundaries, we determined the 
locations of plateaus in the equation of state obtained 
on expansion from a high density columnar crystal 
state. Such plateaus indicate first-order phase transi- 
tions. Phase identification was based on structural prop- 
erties such as the pair-correlation function and nematic 
order parameter, as discussed below. For example, the 
equation of state for L/D = 1 and various values of 
/3i?bond is shown in Figure [31 with the approximate loca- 
tion of first-order phase transitions indicated by horizon- 
tal solid lines. 

Phase diagrams of the sticky-cylinder system for as- 
pect ratios L/D = 0.5, 1 and 2 are shown in FigurePO All 
three phase diagrams share the same key features. First, 
both nematic (N) and columnar LC (C) phases are ob- 
served for large /3£"bond- The range of nematic stability 
(in terms of pressure or packing fraction) increases with 
increasing binding energy. Second, an isotropic-nematic- 
columnar triple point is observed at intermediate bind- 
ing energies. For binding energies below the triple point, 
there is a direct transition from the isotropic phase to 
the columnar LC phase. 

We also find a columnar crystal (X) phase at high 
density, and the L/D = 1 system exhibits a cubatic- 
like phase for small /3-Ebond- The cubatic-like phase was 
previously described by Blaak et al. 80 in the L/D = 
0.9 hard cylinder system. As the focus of this paper is 
on LC phases, we have not investigated the crystalline 
or cubatic phases in detail, and the phase boundaries 
involving these phases should be regarded as schematic. 

Comparison of the equation of state obtained on com- 
pression from an isotropic fluid state with that obtained 
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FIG. 4. Phase diagrams of the sticky cylinder system as a function of binding energy and packing fraction (left) and binding 
energy and pressure (right), for aspect ratios L/D — 0.5 (a), L/D = 1 (b), and L/D = 2 (c). A variety of phases are observed, 
including isotropic (I), nematic (N), columnar liquid crystal (C), columnar crystal (X), and cubatic- like (Cubatic) phases. 
Shaded areas correspond to regions of two-phase coexistence. 



on expansion from a high-density columnar crystal state 
reveals considerable hysteresis (Figure [5]). This hys- 
teresis becomes increasingly pronounced with decreasing 
/3-Ebond, and compression from the low-density isotropic 
state for small /3_Ebond typically leads to a high-density 
disordered (jammed) state, likely due to the dominance 
of hard-core interactions in this limit. For large /3-Ebond, 
the I-N transition occurs at relatively low packing frac- 
tions, which facilitates annealing into a well-ordered ne- 
matic state starting from an isotropic state. Because 
we generally observed slow kinetics of formation of an 
ordered state starting from disordered (isotropic) initial 
configuration, we used isotherms obtained on expansion 



(e.g., Figure [3]) to construct the phase diagrams in Fig- 
ure |4] This implies that the phase boundaries shown are 
lower limits (both in density and pressure), and should 
be considered provisional. In future work we will refine 
the phase diagram using free energy calculations. 



A. Order parameter and correlation functions 

We measured nematic order by calculating the tensor 
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FIG. 5. Hysteresis in the equation of state for L/D = 1 with /3-Bbond = 12 (left) and /3-Ebond = 8 (right). The horizontal lines 
indicate the approximate location of first-order phase transitions observed on expansion from a high-density columnar crystal 
state (solid lines) and on compression from a low-density isotropic fluid state (dashed lines). 
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FIG. 6. Nematic order parameter 5* as a function of pressure 
for L/D = 0.5 and /3E hond = 10. 



where a and j3 refer to cartesian components, i indexes 
the sum over monomers, and 5 a p is the Kronecker delta 
function. The nematic order parameter S is the largest 
eigenvalue of the average order tensor (Q a p), and the ne- 
matic director fi is the corresponding eigenvector. Typ- 
ical results for 5* as a function of pressure are shown in 
Figure § for L/D = 0.5 and pE hond = 10. The order 
parameter jumps abruptly from a small value (S ~ 0.05) 
in the isotropic phase to S ~ 0.7 — 0.8 in the nematic 
phase, and exhibits a further abrupt jump to S *=* 0.9 in 
the columnar LC phase. We list values of the nematic or- 
der parameter in the nematic or columnar liquid crystal 
phases at coexistence (Sn or Sc) with the isotropic phase 
in Table |IJ The values of S are generally large (> 0.7), 
indicating a high degree of orientational order in all LC 
phases investigated. 

To distinguish between the various orientationally or- 



dered phases we computed the three-dimensional pair 
distribution function g(r) and examined positional corre- 
lations as a function of pair separations parallel and per- 
pendicular to the nematic director. The pair distribution 
function g(r) is defined as the normalized probability of 
finding a pair of particles with separation r, 



9( r ) 



1 
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pN EE'c 

* U=i jjti 



(27) 



where p = N/V is the average number density of particles 
and <5(r) is the Dirac delta function. 

Figure f7\ shows example plots of pair distribution func- 
tions. The nematic director defines the z direction 
(z = n), so g(0, y, z) is in a plane containing the nematic 
director and g(x,y,0) is in a plane perpendicular to the 
nematic director. For this set of simulations, L/D = 0.5 
and /3-Et>ond = 12. The columnar LC phase shown in 
Figure IWa) is characterized by a well-defined hexagonal 
lattice of columns and by the absence of long-range cor- 
relations in the z positions of cylinders. The nematic LC 
phase shown in Figure mb) exhibits isotropic and rapidly 
decaying positional correlations in the plane perpendicu- 
lar to n. 

We computed the positional correlation length £y along 
the nematic director in both the nematic and columnar 
LC phases from the one-dimensional pair distribution 
function along a line parallel to the nematic director, 
^(0, 0, z). Figure [8] shows typical results for the columnar 
LC and nematic phases for L/D = 0.5 and /3-Ebond = 12. 
We find exponential decay of correlations to a constant in 
both phases (in the nematic phase g(0, 0, z) approaches 
1 for large z, while in the columnar LC phase we esti- 
mated that g(0,0,z) approaches 3.5 for large z). The 
correlation length determined by fitting the exponential 
decay is £y = 1.3 in the nematic phase at (3PD 3 = 1.2 



10 



(a) 




12 25 38 51 64 77 




S*v 




(b) 




6 12 18 24 30 36 1112 




?\ 0- 



-2- 




-2 2 

X 



FIG. 7. Instantaneous configurations (left) and pair distribution functions in planes containing the nematic director z = n 
(g(0,y, z), center) and perpendicular to the nematic director (g(x,y,0), right) for L/D = 0.5 and j3Ebond — 12. (a) Columnar 
LC phase at J3PD 3 = 1.8, obtained by expansion from a high-density columnar crystal. There are no long-range correlations in 
the z positions of cylinders, as shown by g(0, y, z), and a well-defined hexagonal lattice of columns is apparent in g(x, y, 0). (b) 
Nematic LC phase at /3PD 3 =1.1, obtained by compression from a low-density isotropic fluid. Again, there are no long-range 
correlations in the z positions of cylinders, as shown by g(0, y, z). The plot of g(x, y, 0) reveals isotropic and rapidly decaying 
positional correlations in the plane perpendicular to z, indicative of nematic order. 



and £|| = 1.8 in the columnar LC phase at /3PD 3 = 1.4. 
In the same simulations, the mean aggregate lengths are 
44 in the nematic phase and 352 in the columnar LC 
phase. This gives a ratio of aggregate length to correla- 
tion length that is large: t/&\ ~ 33 in the nematic phase 
and £/£ii ~ 200 in the columnar LC phase. (For more dis- 
cussion of aggregate lengths, see Section IV B| below.) In 
other words, we find that the correlation length can differ 
from the mean aggregate length by 1-2 orders of magni- 
tude. This in turn suggests that it may not be feasible 
to determine the mean aggregation number from X-ray 
scattering experiments that measure correlation lengths, 
such as the experiments of references [321 01] . 



B. Aggregate Statistics 

We measured the aggregation number distribution 
function P(n) and its first moment, the mean aggregation 



number (n), in the isotropic, nematic, and columnar LC 
phases. Here n is the aggregation number (the number 
of monomers in an aggregate) . Aggregates are defined by 
proximity: any pair of cylinders with endpoint separation 
smaller than r c belongs to the same aggregate. 

The use of a finite simulation volume imposes a cutoff 
on the aggregate-length distribution. This finite-size ef- 
fect is particularly evident in the limit of large binding 
energies, where the mean aggregate length can exceed 
the maximum linear dimension of the periodic simula- 
tion cell. To ameliorate this effect, the initial columnar 
crystal configurations used for expansion runs were con- 
structed with an oblique orientation of columns with re- 
spect to the edge vectors of the periodic simulation cell, 
so that a given column traverses the periodic box multi- 
ple times before intersecting itself. With this expedient, 
we are able to measure aggregation number distributions 
over a wide range of conditions, although small finite- 
size artifacts associated with aggregates that connect to 
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FIG. 8. One-dimensional pair distribution function g(0, 0, z) 
along a line parallel to the nematic director in the columnar 
LC phase (blue squares) and the nematic phase (red circles) 
for L/D = 0.5 and /3E bond = 12. 



themselves across the periodic boundaries are observed 
in many of the expansion runs (see Figure pi). 

We illustrate the typical aggregation behavior in Fig- 
ure^! which shows the mean aggregation number (n) as a 
function of pressure (left) and aggregation number distri- 
butions P(n) in the isotropic and nematic states (right) 
for L/D — 0.5 and /3-Ebond = 10, measured upon ex- 
pansion from a high-density columnar crystal state. The 
mean aggregation number jumps discontinuously at the I- 
N transition, and increases strongly with increasing pres- 
sure in the nematic phase. This behavior is qualitatively 
consistent with that observed in previous simulation and 
theoretical studies of semifiexible linear aggregates. The 
distribution P(n) has an exponential dependence on n in 
the isotropic phase, but displays biexponential behavior 
in the nematic phase, i.e., fast exponential decay for small 
n crossing over to slow exponential decay for large n. 
Such biexponential behavior has been noted previously 
by Lii and Kindt [65], who argued that the short length 
scale behavior of P(n) comes from a distinct population 
of short aggregates with relatively low orientational or- 
der. Biexponential behavior in the nematic phase is also 
predicted by the analytic theory; see Figure [TTJ In the 
simulations we also observed a biexponential aggregation 
number distribution in the columnar LC phase. 

Mean aggregation numbers in the isotropic phase at 
coexistence {{n)i) and in the nematic or columnar LC 
phase at coexistence with the isotropic phase {(ti)n or 
(n)c) are listed in Tablefj] The mean aggregation number 
at coexistence tends to increase with increasing /S-Ebond, 
although the trend in the nematic or columnar LC phase 
at coexistence isn't strictly monotonic. We note that 
the mean aggregation numbers at coexistence decrease 
with increasing cylinder aspect ratio for a given jSE'bond, 
possibly because the corresponding packing fractions at 
coexistence decrease with increasing aspect ratio. 



C. Persistence Length 

As discussed above, the persistence length of sticky- 
cylinder aggregates is determined by cylinder-cylinder in- 
teractions. We measured the persistence length from the 
bond orientational correlation function, 



Co(s) = (Ri-R 



^i-\-m 



(28) 



where R^ is a unit vector along the ith bond (the vector 
joining the centers of two neighboring cylinders) in an 
aggregate, the contour length s — m(R), and (R) is the 
average bond length. In the isotropic phase, Co{s) de- 
cays exponentially with increasing m, with a decay con- 
stant that can be identified with the persistence length 
l p , Co(s) = cxp(— s/lp). By fitting this exponential de- 
cay to the simulation data, we determined the persistence 
length. For a given L/D and /3-Ebond, we find that the cal- 
culated persistence length is nearly independent of pres- 
sure within the isotropic phase, which gives us confidence 
that this procedure measures the intrinsic bending rigid- 
ity of aggregates. Table [I] lists the persistence length and 
average bond length as a function of L/D and /SE^ond, 
where both quantities are averaged over isotropic state 
points. 

The persistence length increases approximately lin- 
early with increasing binding energy, and is also approx- 
imately proportional to cylinder aspect ratio L/D for a 
given binding energy. Because of this, the reduced persis- 
tence length l p = lp/ L is reasonably well approximated 
as a universal linear function of /3Eb on d, independent 
of cylinder aspect ratio, as shown in Figure [lOj The 



solid line in Figure 10 is the best linear fit to the data, 
lp = 5.07 + 2.14/32?bond- The ratio of persistence length 
to monomer length approaches a nonzero constant in the 
limit of zero binding energy. This indicates that there 
is a significant entropic contribution to the persistence 
length due to hard-core interactions between neighbor- 
ing cylinders in an aggregate. 

We expect that the precise way that l p varies with 
/3-Ebond to depend strongly on the form of the inter- 
molecular potential. This dependence can be estimated 
by simple dimensional analysis arguments. When two 
bound cylinders interact, thermal fluctuations will typ- 
ically cause their separation to vary to about 1 fc^T 
above their minimum interaction energy E m j n = — Ebond- 
Therefore a typical separation r t = Rj V/^bond- This 
suggests a characteristic angle defined by the distance r t 
and the cylinder radius R: tan 9 = r t /R = \/\f]3E. The 
typical angle is therefore cos 6 = l/-y/l + l/(/3Eb on d)- 
This gives an estimate of the persistence length, based on 



the relation (cos( 



-L/l p 



and assuming /32?b nd ^ L 



of l p — 2/3i?bond- In other words, the persistence length 
should vary linearly with /3 Ebond, with a slope of 2, as 
observed in simulations. 

For other forms of the potential, this dependence 
would change. Using the same dimensional analysis ar- 
gument, we would predict that for an exponent of 7 in 
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FIG. 9. Mean aggregation number (n) as a function of pressure (left) and aggregation number distributions P(n) (right) for 
L/D — 0.5 and /3-Ebond = 10, measured on expansion from a high-density columnar crystal state. The nematic phase P(n) 
displays fast exponential decay for small n crossing over to slow exponential decay for large n (bi-exponential distribution). 
The small peak near n = 20 in the nematic aggregation number distribution is a finite-size artifact, associated with aggregates 
that connect to themselves across the periodic boundaries. 



TABLE I. Properties of the sticky-cylinder system measured in NPT MC simulations. Columns contain the cylinder aspect 
ratio L/D, the dimensionless binding energy /3_Ebond, the number of particles N, the measured persistence length l p , the mean 
intra-aggregate bond length (R), the mean aggregation number in the isotropic phase at coexistence (n)j, the mean aggregation 
number in the nematic or columnar LC phase at coexistence with the isotropic phase ((n)jv or (n)c), and the nematic order 
parameter in the nematic or columnar phase at coexistence with the isotropic phase (Sn or Sc)- 
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FIG. 10. Ratio of persistence length to monomer length 
l p = lp/L as a function of /9-Bbond for aspect ratios L/D = 
0.5 (circles), L/D = 1 (squares), and L/D — 2 (triangles). 
The data for all three aspect ratios is well approximated by 
a universal linear function, l p — 5.07 + 2.14/3_Bbond, which 
represents the best linear fit to the data (solid line). 



the interaction potential (Eqn. (fl])), the scaling is l p = 
2(/3i?bond) 2 ^ 7 - In the limit of a square well (7 — » 00) the 
persistence length would become independent of /3£t>ond- 
Because the phase behavior of aggregates is strongly de- 
pendent on the persistence length, this implies a strong 
dependence of the phase behavior on the form of the 
monomer-monomer interaction potential, in addition to 
the obvious dependence on the binding energy. In the 
future it would be interesting to investigate the interplay 
of these two effects, focusing how sticky-cylinder phase 
behavior varies with the form of the interaction potential. 
The aggregates in this study are relatively flexible: the 
ratio of persistence length to mean aggregate length for 
the nematic phase at coexistence ranges from w 0.8 for 
L/D = 0.5 and /3E hond = 12 to « 2.9 for L/D = 2 and 

/3-Ebond = 12. 



V. THEORY-SIMULATION COMPARISON 

We compared our analytic theory results both to the 
simulations of Lii and Kindt 65 and to our simulations. 
These comparison illustrate the key role of aggregate flex- 
ibility in determining the theory-simulation agreement. 

First, we calculated the aggregate length distributions 
and order parameter as a function of length to compare 



with the simulation results of Lii and Kindt (Figure 11 1 



Lii and Kindt considered two persistence length values, 
lp = 1000 and l p = 100 65J . Note that l p = l p /L is the 
persistence length divided by the monomer length. (In 
Lii and Kindt's work, the monomers are spheres.) For 
this figure, we used the the association constant K — 
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5000. 



A significant and unexpected result is the calculated bi- 
exponential aggregate length distribution in the nematic 
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FIG. 11. Theoretical predictions of aggregate length and 
order parameter distributions. Top, aggregate length distri- 
butions at I-N coexistence of the isotropic (I, dashed curves) 
and nematic (N, solid curves) phases with l p — 1000 and 
l p = 100 (subscripts). Note that in this figure the predictions 
of the simple second-virial approximation and the Parsons- 
Lee approximation overlap. Bottom, length-dependent ag- 
gregate order parameter at I-N coexistence in the nematic 
phase with l v — 1000 and l v — 100 (subscripts). Solid curves: 
simple second-virial approximation. Dashed curves: Parsons- 
Lee approximation. In both cases, the association constant 
K = e^ Bbond = 5000, for comparison with the simulations of 
Lii and Kindt 1651. 



phase. This distribution arises naturally from the free- 
energy minimization (Figure 



11 



top). To our knowledge, 
our work is the first to find this biexponential distribu- 
tion in analytic theory. Lii and Kindt's initial analytic 
work found only a single exponential distribution |65j : 
in later work they found that assuming a biexponential 
distribution improved the theory-simulation agreement 

Overall we find good quantitative agreement between 
the predictions of the analytic theory for aggregate length 
distributions and order parameter as a function of aggre- 
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FIG. 12. Theoretical predictions of mean aggregation num- 
ber in the isotropic and nematic phases at I-N coexistence. 
(A), l p = 1000. (B), l p = 100. Comparison of theory (solid 
and dashed curves) to simulation results of Lii and Kindt 
|65| (points). Curves denote borders of regions of I-N co- 
existence. Solid curves: simple second-virial approximation. 
Dashed curves: Parsons-Lee approximation. Points and bars: 
simulation results of Lii and Kindt |65| . Note that for l p = 100 
(B) Lii and Kindt gave ranges of measured values, which we 
plot using horizontal bars, while for l p = 1000 (A) Lii and 
Kindt gave only mean values, which we plot using points. 



gate length (compare Figure 11 to Figure 2 of reference 
[55]). The aggregate length distribution is exponential 
in the isotropic phase and biexponential in the nematic 
phase. The nematic phase also shows significantly longer 
aggregates. Stiffer aggregates also tend to be longer 
in the nematic phase. The average alignment depends 
strongly on aggregate length in the nematic phase, with 
short aggregates only weakly aligned and a rapid increase 
in the length-dependent order parameter with aggregate 
length. 

Next we compared our predictions for the aggregate 
length and nematic order parameter at FN coexistence to 
the simulations of Lii and Kindt for a range of monomer 
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FIG. 13. Comparison of nematic order parameter at I-N 
coexistence predicted by theory (dashed curves) to simulation 
results of Lii and Kindt 65 (points). Persistence lengths 
l p = 1000 and l p = 100 indicated by subscripts. 



binding energies -©bond- In Figure [12] we plot the av- 
erage aggregation number in the isotropic and nematic 
phases at FN coexistence for a range of values of i^bond- 
The curves are predictions of our analytic theory and 
the bars and points are values from the simulations of 
Lii and Kindt [55J. Note that for the shorter persistence 
length l p = 100 Lii and Kindt gave ranges of measured 
values, which we plot using horizontal bars, while for 
the larger persistence length l p = 1000 Lii and Kindt 
gave only mean values, which we plot using points. (The 
same holds for the two figures below in which we com- 
pare to Lii and Kindt's work). This result illustrates the 
strong enhancement of aggregation in the nematic phase. 
The theory clearly captures the correct qualitative trends 
of the simulations: as the binding energy increases, the 
aggregates become longer in both the isotropic and the 
nematic phases. Quantitatively agreement is reasonable 
for the larger persistence length, but for the shorter per- 
sistence length l p = 100 the predicted values of (n) are 
off by approximately 50%. This illustrates the decreased 
quantitative accuracy of the theory as the aggregates be- 
come less flexible. 

The order parameter in the nematic phase at FN co- 
existence is shown in Figure |13[ Here our calculations 
for the order parameter agree well with the simulation 
results. As expected, the stiffer, longer aggregates show 
a higher average order parameter that is less sensitive to 
the binding energy. 

Our final comparison with the simulation results of Lii 
and Kindt is the phase diagrams shown in Figure |14[ 
The analytic theory agrees well with the l p — 1000 simu- 
lation results and less well with the l p — 100 simulation 
results. For l p = 100 the analytic theory gives errors 
of up to about 50% in the predicted volume fraction <p 
of the coexistence region. Perhaps surprisingly, we find 
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FIG. 14. Comparison of theory (solid and dashed curves) 
to simulation results of Lii and Kindt ^65 (points). (A), 
l p — 1000. (B), l p = 100. Curves denote borders of re- 
gions of I-N coexistence. Solid curves: simple second-virial 
approximation. Dashed curves: Parsons-Lee approximation. 



that the simple second-virial approximation (solid curves 
in Figure [m] ) are closer to the simulation results than 
the results that include the Parsons-Lee approximation 
(dashed curves). 

In our simulations, the ratio of persistence length to 
monomer length l p = L/D is significantly less than in 
the work of Lii and Kindt; we find values l p « 20 — 30 
in the I-N coexistence range (Table II]). The ratio of mean 
aggregate length to persistence length £/l p w 0.1 — 1 in 
this range (Table|l|. Because our analytic theory assumes 
nearly rigid aggregates (£/l p <C 1) we expect that the the- 
ory will not show good quantitative agreement with the 
simulations in this regime. We note that in our simu- 
lations the persistence length is not a control parame- 
ter but is determined by the balance of binding energy 
and entropy of the aggregates; we found empirically that 
lp = 5.07 + 2.14:f3E hond (figure 10). We used this empiri- 
cal relationship in the analytic theory. 




FIG. 15. Comparison of I-N coexistence in theory (curves) 
and simulation (points). Top, aggregation number. Botton, 
phase diagarm Curves denote borders of regions of I-N co- 
existence; gray: onset of I-N coexistence; black: upper limit 
of I-N coexistence. Solid curve: L/D — 2. Dashed curves: 
L/D — 1. Dotted curves: L/D = 0.5. In doing the calcu- 
lations we used the empiricial persistence length found from 
simulations l p = 5.07 + 2.14/3i?bond. Triangles: L/D — 2. 
Squares: L/D = 1. Circles: L/D = 0.5. 



In Figure 15 we show the aggregation number and vol- 
ume fraction at the borders of I-N coexistence. While 
the qualitative trends in the simulations arc correctly 
captured by theory, there is significant quantitative dis- 
agreement. The analytic theory predicts significantly 
larger aggregation numbers and significantly lower vol- 
ume fractions than found in the simulations. The biggest 
discrepancy is with the aggregation numbers, where the 
analytic theory predictions are larger than found in sim- 
ulation by factors of 2-5. For the volume fraction, the 
discrepancy is a factor of 1/3-4. However, we do note 
that although the position of the phase boundaries is not 
correctly predicted by our theory, their slope is. We note 
that the slope is affected by the variation of the persis- 
tence length l p with /3-Ebond : including this variation in 
the analytic theory is necessary to reproduce the slope 
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FIG. 16. Comparison of order parameter in theory (curves) 
and simulation (points) in the nematic phase at I-N coexis- 
tence. Solid curve: L/D — 2. Dashed curves: L/D — 1. 
Dotted curves: L/D = 0.5. In doing the calculations we 
used the empiricial persistence length found from simulations 
l v = 5.07 + 2.U/3E hond . Triangles: L/D = 2. Squares: 
L/D = 1. Circles: L/D = 0.5. 



seen in the simulations. 

We also find significant quantitative discrepancy in the 
order parameter in the nematic phase at coexistence (Fig- 
ure 



16 1 . The theory predicts lower values of S than found 



in the simulations by almost a factor of 2. 



VI. CONCLUSION 

Recent experiments on the coupled aggregation and 
liquid-crystal ordering of chromonic liquid crystals and 
short double-stranded nucleic acids motivated us to con- 
sider coarse-grained simulation and theory of this prob- 
lem. In nucleic acid systems in particular, the monomer 
aspect ratio can be varied by varying the number of 
base pairs per monomer; however, the effects of varying 
monomer aspect ratio on the aggregation/LC ordering 
problem have not been extensively studied in theoretical 
work. 

In this paper we used both Monte Carlo simulation and 
analytic theory to study aggregation and liquid-crystal 
ordering of a simple model of sticky cylinders. The model 
assumes hard cylinders (of length L and diameter D, so 
the aspect ratio of the monomers can easily be varied) 
that experience attractive interactions when they stack 
end to end. 

We determined approximate phase diagrams by using 
calculations of the equation of state in expansion and 
compression runs; we then studied the aggregate length 
distributions, order parameter, and correlation functions 
in the different phases. We find isotropic, nematic, 



columnar LC, columnar crystal, and cubatic-like phases 
are possible for this system. This family of phases is simi- 
lar to the results of Hentschke and colleagues [751 EZ] but 
does show differences with the work of Lii and Kindt, 
who only observed isotropic and nematic phases [65] . 
Similarly to previous work we observe the disappearance 
of the nematic phase (as the monomer binding energy 
drops or temperature increases) at an isotropic-nematic- 
columnar triple point g2 S3 [Ml E3 [77] ■ These phase 
diagrams bear a strong qualitative resemblance to those 
of chromonic liquid crystals, which show a strong de- 
crease in nematic phase width with increasing temper- 
ature 02 gig. The location of the I-N-C triple point 
appears to be sensitive to the monomer aspect ratio, oc- 
curring at lower /3E bond for larger L/D. 

In the simulations we studied three aspect ratios: 
L/D — 0.5, 1, and 2. The same qualitative phase behav- 
ior is present for all aspect ratios. The key differences 
are that, first, we only observe the cubatic-like phase for 
L/D = 1. This is consistent with previous work on hard 
cylinders |80j . Second, we do not observe the columnar 
crystal phase for the lowest aspect ratio of 0.5. Naively 
this might be understood because having columns made 
of a larger number of shorter monomers makes the in- 
crease in entropy allowed by disorder along the columns 
more favorable. Third, as noted above, the nematic phase 
persists to lower /3-Ebond when the aspect ratio is larger. 
Qualitatively, this makes sense. We find that the mean 
aggregate length at I-N coexistence increases with in- 
creasing monomer aspect ratio, which suggests (as in On- 
sager theory) that the density at coexistence should also 
decrease with increasing monomer aspect ratio. We do 
observe this decrease, although the trend is not extremely 
strong. We also find that the persistence length increases 
linearly with increasing monomer aspect ratio; therefore 
aggregate stiffness increases with increasing aspect ra- 
tio (assuming fixed binding energy) . Both increasing ag- 
gregate length and increasing aggregate stiffness tend to 
favor nematic order, and therefore tend to broaden the 
nematic region. 

Our results give evidence that the phase diagram de- 
pends on the the form of the interaction potential, not 
just the binding energy between monomers. The phase 
behavior is sensitive to the flexibility of the aggregates, 
and the aggregate persistence length varies with temper- 
ature in a way that is controlled by the interaction po- 
tential. This suggests that the temperature dependence 
of the I-N phase boundary as well as the location of the 
I-N-C triple point arc influenced by the implicit depen- 
dence of l p on /3i?bond- This is an interesting direction 
for future work. 

Our simulation results demonstrate a remarkable dif- 
ference between the correlation length along aggregates 
and the length of aggregates in the nematic and colum- 
nar LC phases. Indeed, the aggregate length is typ- 
ically 10-100 times larger than the density correlation 
length along the columns. This has important implica- 
tions for attempts to measure aggregate length distri- 
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butions: X-ray scattering experiments, which measure 
correlation lengths, may not be able to probe aggregate 
lengths Hi]. 

In our computations the aggregate persistence length 
is determined from the interaction potential of the 
monomers; we found a universal (independent of L/D) 
linear dependence of l p = l p /L on £t,ond- 

Our extensions of previous theoretical work suggest 
that quantitative agreement between analytic theory and 
simulation for the I-N transition of aggregates is within 
reach. Starting with the relatively simple second-virial 
model of aggregation and nematic order of van der Schoot 
and Cates [S3], we added (1) the Parsons-Lee approxi- 
mation for higher-density systems, and (2) expansion of 
the angular trial function in Legendre polynomials and 
numerical solution of the resulting algebraic equations. 
This combination significantly improves the quantitative 
agreement between the analytic theory and the simula- 
tions of Lii and Kindt [65 . 

Our analytic theory predicts the biexponential distri- 
bution of aggregate lengths in the nematic phase. Previ- 
ous analytic work of Lii and Kindt assumed a biexponen- 
tial distribution in the nematic phase, and based on this 
assumption found good agreement with their simulation 
results [66 . Here the biexponential distribution arises 
directly from the free-energy minimization and does not 
have to be added to the theory as an assumption. 

The comparison between simulation and theory high- 
lights the role of aggregate flexibility in controling the 
phase diagram. In our analytic theory, the aggregates 
are treated as nearly rigid. Therefore, our theoretical 
results match the simulation results best where the ag- 
gregate persistence length is long (l p — 1000). In this 
regime, the quantitative agreement between theory and 
simulation is remarkably good, particularly in the pre- 
diction of the phase boundaries, the order parameter, 
and the biexponential aggregate length distribution in 
the nematic phase. When comparing to the simulations 
of Lii and Kindt with shorter aggregate persistence length 
(l p = 100) or our simulations (l p sw 20 — 30), the quanti- 
tative agreement between analytic theory and simulation 
decreases. Future improvement to the theory will require 
improved treatment of aggregate flexibility. 

A major direction for future research is the compari- 
son of our results to experimental results on chromonic 
and nucleic-acid liquid crystals, and refinement of both 
the simulation model and the analytic theory based on 
the comparison. In the future, improvements to the sim- 
ulation model and the analytic theory may allow true 
quantitative agreement between experiment, simulation, 
and theory. 



Appendix A: Free energy minimization in the 
nematic phase 

Here we outline the free energy minimization and nu- 
merical solution for the function M(x) introduced in 



Eq. (18 1. Plugging the trial function Eq. (18) into the 
normalization condition Eq. (f5j) we arrive at the normal- 
ization condition for M(x) = M(x)/M 



dxM 2 (x) = 2. 



(Al) 



Minimizing the free energy /("' (Eq. (19)) with respect 
to M(x) subject to the normalization constraint gives 



= 0, 



S(f (n) +Xj\dxM 2 (x)) 
SM(x) 

leading to the equation for M(x), 



M{x) = \M 2 {x) - ^M(x)d 2 x M{x) 

dip 



(A2) 



4r0M o 



dx'M 2 (x)M 2 (x')K(x,x'). (A3) 



Integrating Eq. ( A3 ) over x and using the normalization 



condition we obtain the equation for the Lagrange mul- 
tiplier A: 



A 



1 



Mo 

Of i> 



h 



4r</>M 



I, 



where the three integrals are 
h = / dxM(x), 



dxM(x)Q%M(x), 



(A4) 

(A5) 
(A6) 



-l 
i 



[...= I dx 1 dx 2 K(x 1 ,x 2 )M 2 {x 1 )M 2 {x 2 ). (A7) 



To solve equation ( A3 ) and determine the function 



M(x) that minimizes the free energy we seek a series so- 
lution: expansion of M(x) in Legendre polynomials turns 
equation ( A3 1 into a system of coupled algebraic equa- 



tions that can straightforwardly be truncated and solved 
numerically. We expand M(x) in even Legendre polyno- 
mials P 2k (x) (note that the nematic phase is symmetric 
under x — > —x, so only even Legendre polynomials are 
allowed): 



M{x) = ^m 2k P 2k {x) 



(A8) 



k=0 



Plugging this series expansion into Eq ( A3 1 gives the 



equations for the coefficients m 2k , after using the orthog- 
onality condition and completeness relation for Legendre 
polynomials: 



m 2k 



4k + 1 



i "2fc + ^2 { ( 4i + !) fc 2j Ai a 2l a 2j 
(2i)(2i + l)A 2 m2im2j\l2k,2i,2j ■ (A9) 
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The dimensionless parameters Ai and A2 are 



Ai 



A, 



Ar<j)M Q 



2M n 



o / . 



(A10) 
(All) 



The coefficients k 2 j come from the expansion of the in- 
teraction kernel 



2t, 



. #1 d(j> 2 , 
K(x 1 ,x 2 ) = I --^_7^_7\ sm "f\ 



2tt 2tt 



(A12) 



in Legendre polynomials, 



K{xi,X 2 ) = ^k 2n P 2n (x!)P 2n (x 2 ), (A13) 



n=0 



with fco — 7r /4, k 2 — — 57r/32 and 

7r(4n + l)(2n-3)!!(2n-l)!! 

2n ~ 2 2 «+ 2 n!(n+l)! 

The coefficients a 2n are 

a 2n = 2 2_^ TO2a ™2b ^2a,2b,2m 



(n> 1). 
(A14) 

(A15) 



other cases 

l2k,2i,2j = 

(2k + 2i- 2j)\(2k -2i + 2j)\(-2k + 2i + 2j)\ 



(2k + 2i + 2j + 1)\ 
(k + i + j)) 



(k + i- j)\(k - i + ])\(-k + i + j)\ 



(A17) 



The three integrals in Eqs. (A5)-(A7) become, after 
plugging in the series expansion Eq. (A8), 



h = 2ra , 

'* = -£ 

h 



4n(2n+l), 2 



4n + 1 

n 

/ J k 2n (a2n) 2 - 



(A18) 
(A19) 

(A20) 



We truncated and solved the system of algebraic equa- 
tions dA4h, (|A9| and (|Al8|)-(|A20b numerically in Mat- 



lab. We typically kept 25-30 expansion coefficients for 
the function M(x), depending on how well the numerical 
solution converged. 

In terms of the expansion coeficients m 2 k the mean 
aggregation number (n) in the nematic phase is 



2Mn 



A/n 



J_ 1 dxM(x) m 
and the order parameter is 
j\dxP 2 (x)M 2 (x 



(A21) 



S = 



j\dxM 2 (x) 



= y^ j m2iTn2jh,2i,2j- (A22) 



and the integral of 3 Legendre polynomials is 



I 



2k,2i,2j 



1 f 1 

- / dxP 2k (x)P 2l (x)P 2] (x), (A16) 



This integral is known exactly [S5]: it is equal to when 
k + i — j < or k — i + j < or — k + i + j < 0. In all 
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